Dynamical instability of new-born neutron stars as sources of gravitational radiation 
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The dynamical instability of new-born neutron stars is studied by evolving the linearized hy- 
drodynamical equations. The neutron stars considered in this paper are those produced by the 
accretion induced collapse of rigidly rotating white dwarfs. A dynamical bar-mode (m = 2) in- 
stability is observed when the ratio of rotational kinetic energy to gravitational potential energy 
[3 of the neutron star is greater than the critical value fid « 0.25. This bar-mode instability leads 
to the emission of gravitational radiation that could be detected by gravitational wave detectors. 
However, these sources are unlikely to be detected by LIGO II interferometers if the event rate is 
less than 10 -6 per year per galaxy. Nevertheless, if a significant fraction of the pre-supernova cores 
are rapidly rotating, there would be a substantial number of neutron stars produced by the core 
collapse undergoing bar-mode instability. This would greatly increase the chance of detecting the 
gravitational radiation. 



PACS: 04.30.Db, 95.30.Sf, 97.60.Jd 

I. INTRODUCTION 

Neutron stars are believed to form from the core col- 
lapse of massive stars and the accretion induced collapse 
of massive white dwarfs. If the stellar core or white dwarf 
is rotating, conservation of angular momentum implies 
that the resulting neutron star must rotate very rapidly. 
It has been suggested Q that such a rapidly rotating star 
may develop a non-axisymmetric dynamical instability, 
emitting a substantial amount of gravitational radiation 
which might be detectable by gravitational wave obser- 
vatories such as LIGO, VIRGO, GEO and TAMA. 

Rotational instabilities arise from non-axisymmetric 
perturbations having angular dependence e lmip , where (p 
is the azimuthal angle. The m = 2 mode is called the 
bar mode, which is usually the strongest mode for stars 
undergoing instabilities. There are two types of instabili- 
ties. A dynamical instability is driven by hydrodynamics 
and gravity, and it develops on a dynamical timescale, i.e. 
the timescale for a sound wave to travel across the star. 
A secular instability, on the other hand, is driven by vis- 
cosity or gravitational radiation reaction, and its growth 
time is determined by the relevant dissipative timescale. 
These secular timescales are usually much longer than 
the dynamical timescale of the system. 

In this paper, we focus on the dynamical instabili- 
ties resulting from the new-born neutron stars formed 
from accretion induced collapse (AIC) of white dwarfs. 
These instabilities occur only for rapidly rotating stars. 
A useful parameter to characterize the rotation of a star 
is (3 = T/\W\, where T and W are the rotational ki- 
netic energy and gravitational potential energy respec- 
tively. It is well-known that there is a critical value 
(3d so that a star will be dynamically unstable if its 
[3 > [3d- For a uniform density and rigidly rotating 
star, the Maclaurin spheroid, the critical value is deter- 
mined to be (3d ~ 0.27 Numerous numerical sim- 
ulations using Newtonian gravity show that (3d remains 



roughly the same for differentially rotating polytropes 
having the same specific angular momentum distribution 
as the Maclaurin spheroids P-p"l| . However, (3d can take 
values between 0.14 to 0.27 for other angular momentum 
distributions p2|,p 13 (the lower limit (3d = 0.14 is ob- 
served only for a star having a toroidal density distribu- 
tion, i.e. the maximum density occurs off the center p3|). 
Numerical simulations using full general relativity and 
post-Newtonian approximations suggest that relativistic 
corrections to Newtonian gravity cause (3d to decrease 
slightly 

Most of the stability analyses to date have been car- 
ried out by assuming that the star rotates with an ad 
hoc rotation law or using simplified equations of state. 
The results of these analyses might not be applicable to 
the new-born neutron stars resulting from AIC. Recently, 
Fryer, Holz and Hughes |l7|] carried out an AIC simula- 
tion using a realistic rotation law and a realistic equation 
of state. Their pre-collapse white dwarf has an angular 
momentum J = 10 49 g cm 2 s _1 . After the collapse, the 
neutron star has (3 less than 0.06, which is too small for 
the star to be dynamically unstable. However, they point 
out that if the pre-collapse white dwarf spins faster, the 
resulting neutron star could have high enough (3 to trigger 
a dynamical instability. They also point out that a pre- 
collapse white dwarf could easily be spun up to rapid ro- 
tation by accretion. The spin of an accreting white dwarf 
before collapse depends on its initial mass, its magnetic 
field strength and the accretion rate, etc. jig] . 

Liu and Lindblom jl9| (hereafter Paper I) in a recent 
paper construct equilibrium models of new-born neutron 
stars resulting from AIC based on conservation of spe- 
cific angular momentum. Their results show that if the 
pre-collapse white dwarfs are rapidly rotating, the result- 
ing neutron stars could have (3 as large as 0.26, which is 
slightly smaller than the critical value (3d for Maclaurin 
spheroids. However, the specific angular momentum dis- 
tributions of those neutron stars are very different from 
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that of Maclaurin spheroids. So there is no reason to be- 
lieve that the traditional value (3d = 0.27 can be applied 
to those models. 

The purpose of this paper is first to determine the crit- 
ical value (3d for the new-born neutron stars resulting 
from AIC, and then estimate the signal to noise ratio 
and detectability of the gravitational waves emitted as 
a result of the instability. We do not intend to provide 
an accurate number for the signal to noise ratio, which 
requires a detailed non-linear evolution of the dynamical 
instability. Instead, we use Newtonian gravitation the- 
ory to compute the structure of new-born neutron stars. 
Then we evolve the linearized Newtonian hydrodynami- 
cal equations to study the star's stability and determine 
the critical value (3d- Relativistic effects are expected to 
give a correction of order (v/c) 2 , which is about 8% for 
the rapidly rotating neutron stars studied in this paper. 
Here v is a typical sound speed inside the star and c is 
the speed of light. 

This paper is organized as follows. In Sec. ||, we apply 
the method described in Paper I to construct a number 
of equilibrium neutron star models with different values 
of (3. In Sec. Ill, we study the stability of these mod- 



els by adding small density and velocity perturbations 
to the equilibrium models. Then we evolve the pertur- 
bations by solving linearized hydrodynamical equations 
proposed by Toman et al [ p9| . From the simulations, 
we can find out whether the star is stable, and determine 
the critical value (3d- In Sec. we estimate the strength 
and signal to noise ratio of the gravitational waves emit- 
ted by this instability. In Sec. [V| we estimate the effects 
of a magnetic field on the stability results. Finally, we 
summarize and discuss our results in Sec. [v|. 



II. EQUILIBRIUM MODELS 

In this section, we describe briefly how we construct 
new-born neutron star models from the pre-collapse 
white dwarfs. A more detailed description is given in 
Paper I. 



A. Pre-collapse white dwarf models 

We consider two types of pre-collapse white dwarfs: 
those made of carbon-oxygen (C-O) and those made of 
oxygen- neon-magnesium (O-Ne-Mg). The collapse of a 
massive C-0 white dwarf is triggered by the explosive 
carbon burning near the center of the star p0| , pl |. The 
central density of the pre-collapse C-0 white dwarf must 
be in the range 6 x 10 9 g cm~ 3 < p c < 10 10 g cm~ 3 in 
order for the collapse to result in a neutron star, rather 
than exploding as a Type la supernova [^2| . The collapse 
of a massive O-Ne-Mg white dwarf, on the other hand, 
is triggered by electron captures by 24 Mg and 20 Ne when 
the central density reaches 4 x 10 9 g cm~ 3 pg,Bl|. 



We construct three sequences of pre-collapse white 
dwarfs, with models in each sequence having different 
amounts of rotation. Sequences I and II correspond to 
C-0 white dwarfs with central densities p c — 10 10 g cm~ 3 
and p c = 6 x 10 9 g cm -3 respectively. Sequence III is for 



O-Ne-Mg white dwarfs with p c = 4 X 10 g cm 
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All 



white dwarfs are assumed to rotate rigidly, because the 
timescale for a magnetic field to suppress differential ro- 
tation is much shorter than the accretion timescale (see 
Sec.|v3). 

The pre-collapse white dwarfs constructed in this sec- 
tion are described by the equation of state (EOS) of a 
zero-temperature ideal degenerate electron gas with elec- 
trostatic corrections derived by Salpeter |23|] . At high 
density, the pressure is dominated by the ideal degener- 
ate Fermi gas with electron fraction ZJA = 0.5 that is 
suitable for both C-0 and O-Ne-Mg white dwarfs. Elec- 
trostatic corrections, which depend on the white dwarf 
composition through the atomic number Z, contribute 
only a few percent to the EOS for the high density white 
dwarfs considered here. 

Equilibrium models are computed by Hachisu's self- 
consistent field method |2j] , which is an iteration scheme 
based on the integrated Euler equation for hydrostatic 
equilibrium: 



h + <t>- — £T=C, 



(1) 



where fi is the rotational angular frequency of the star; 
C is a constant; w is the radius from the rotation axis; 
h is the specific enthalpy, which is related to the density 
p and pressure P by 



Jo P 



(2) 



The gravitational potential <& satisfies the Poisson equa- 
tion 



V 2 $ = AnGp 



(3) 



where G is the gravitational constant. The self-consistent 
field method determines the structure of the star for fixed 
values of two adjustable parameters. In Ref. [Q, the 
maximum density and axis ratio (the ratio of polar to 
equatorial radii) are the chosen parameters. However, it 
is more convenient to choose the central density p c and 
equatorial radius R e as the two parameters for the mod- 
els studied here. 
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TABLE I. Properties of pre-collapse white dwarfs. Here S7 
is the rotational angular frequency; £l m is the maximum ro- 
tational angular frequency of the white dwarf in the sequence 
without mass-shedding; R e , R p , M, J and (3 are respectively 
the equatorial radius, polar radius, mass, angular momentum 
and the ratio of rotational kinetic to gravitational potential 
energies. 

Sequence I: C-0 white dwarfs with p c = 10 10 g cm" 3 





Re 


Re/ Rp 


M/Mq 


J 


13 






(km) 






(g cm 2 s" 1 ) 






0.20 


1310 


1.01 


1.40 


5.14 x 10 48 


5.36 x 10" 


-4 


0.65 


1400 


1.09 


1.42 


1.81 x 10 49 


6.08 x 10" 


-3 


0.84 


1500 


1.19 


1.44 


2.43 x 10 49 


1.03 x 10" 


-2 


0.93 


1600 


1.27 


1.46 


2.80 x 10 49 


1.30 x 10" 


-2 


1.00 


1895 


1.52 


1.47 


3.12 x 10 49 


1.55 x 10" 


-2 



Sequence II: C-0 white dwarfs with 
p c = 6 x 10 10 g cm" 3 





Re 


Re/ Rp 


M/Mq 


J 


P 






(km) 






(g cm 2 s" 1 ) 






0.23 


1517 


1.01 


1.39 


3.31 x 10 48 


7.73 x 10" 


-4 


0.64 


1610 


1.09 


1.42 


1.95 x 10 49 


6.11 x 10" 


-3 


0.84 


1740 


1.19 


1.44 


2.75 x 10 49 


1.12 x 10" 


-2 


0.93 


1847 


1.28 


1.45 


3.13 x 10 49 


1.39 x 10" 


-2 


1.00 


2189 


1.52 


1.46 


3.51 x 10 49 


1.66 x 10" 


-2 



Sequence III: O-Ne-Mg white dwarfs with 
p c = 4 x 10 10 g cm" 3 





Re 


Re/ Rp 


M/Mq 


J 


P 






(km) 






(g cm 2 s" 1 ) 






0.23 


1692 


1.01 


1.38 


7.01 x 10 48 


7.90 x 10" 


-4 


0.62 


1791 


1.09 


1.40 


2.05 x 10 49 


6.23 x 10" 


-3 


0.86 


1956 


1.20 


1.42 


3.03 x 10 49 


1.23 x 10" 


-2 


0.96 


2156 


1.34 


1.44 


3.59 x 10 49 


1.62 x 10" 


-2 


1.00 


2441 


1.52 


1.45 


3.80 x 10 49 


1.77 x 10" 


-2 
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The accuracy of the equilibrium models can be mea- 
sured by the quantity 



2T + W + 311 



w 



(4) 



which should be equal to zero according to the Virial 
theorem (see e.g. [25j). Here T is the rotational ki- 
netic energy; W is the gravitational potential energy and 
II = J Pd 3 x. The values of e are of order 10 -7 for all 
the pre-collapse white dwarf models calculated in this 
section. Table | shows some properties of several pre- 
collapse white dwarfs in the three sequences. Each se- 
quence terminates when the rotational angular frequency 
of the white dwarf reaches a critical value f2 m so that 
the mass-shedding occurs on the equatorial surface of the 
star. The values of Cl m are 5.37 rad s _1 for Sequence I, 
4.32 rad s _1 for Sequence II, and 3.65 rad s _1 for Se- 
quence III. 




0.2 0.4 0.6 O.f 

FIG. 1. The normalized specific angular momentum j as 
a function of the cylindrical mass fraction 7n ro for the white 
dwarf model in sequence III with fi/fl m = 0.964. 

Fig. |l| shows the normalized specific angular momen- 
tum 



j = —£lw 2 
J 

as a function of the cylindrical mass fraction 



(5) 



m ro = — j d—'zu' j (/:'/'( zz '. :' ) ((>) 



for a typical pre-collapse model. Here M and J are the 
total mass and angular momentum of the star respec- 
tively. The specific angular momentum defined in Eq. (|5|) 

is normalized so that J Q j{m^) dm^ = 1. The curves for 
other models are similar. The spike near m ro = 1 is due 
to the high degree of central condensation of the white 
dwarf density (Paper I). 



B. Collapsed objects 

The gravitational collapse of a massive white dwarf is 
halted when the core density reaches nuclear density. The 
core bounces back and settles down into hydrodynamical 
equilibrium in a few milliseconds. A hot (T > 20 Mev) 
and lepton-rich proto-neutron star is formed. After about 
20 s, neutrinos carry away most of the energy and the star 
cools down to a cold neutron star. The hot proto-neutron 
stars are less compact and have f3 smaller than 0.14 (Pa- 
per I). They are thus expected to be dynamically stable. 
Hence in this paper, we focus on the stability of the cold 
neutron stars shortly after the cooling. 

We assume that (1) the neutron stars are axisymmet- 
ric and are in rotational equilibrium with no meridional 
circulation; (2) viscosity can be neglected; (3) no ma- 
terial is ejected during the collapse. Under these as- 
sumptions, it is easy to prove that (see e.g. |l9| , ^5|l ) the 
specific angular momentum j of the collapsed star as a 
function of cylindrical mass fraction m ro is the same as 
the pre-collapse white dwarf. Hence the structure of the 
new-born neutron stars can be constructed by computing 
models with the same masses, angular momentum and 
i(m ro )-distributions as the pre-collapse white dwarfs. 

We adopt the Bethe-Johnson EO S p6[ for densities 
above 10 14 g cm" 3 , and BBP EOS (27j for densities in 
the range 10 11 — 10 14 g cm' 3 . The EOS for densities 
below 10 11 g cm -3 is joined by that of the pre-collapse 
white dwarfs. 

We construct the equilibrium models by the numerical 
method proposed by Smith and Centrella |^] , which is a 
modified version of Hachisu's self-consistent field method 
so that j(?n ro ) can be specified. The iteration scheme is 
based on the integrated Euler equation ([]]) written in the 
form 



dvo 



, .finite') 



(7) 



where J and M are the total angular momentum and 
mass of the star respectively. As before, two parameters 
have to be fixed in the iteration procedure. We choose to 
fix the central density p c and equatorial radius R e . Since 
the correct p c and R e are not known beforehand, we have 
to vary these two quantities until the equilibrium model 
has the same J and M as the pre-collapse white dwarf. 

The standard iteration algorithm described in 
Refs. |Q and [M fails to converge when the star be- 
comes very flattened. This problem is fixed by a modi- 
fied scheme proposed by Pickett, Durisen and Davis 
in which only a fraction of the revised density (or en- 
thalpy) pi +1 , i.e. p' i+l = (1 - Qpi+i + (Pi, is used for 
the next iteration. Here C < 1 is a parameter controlling 
the change of density. A value of C > 0.95 has to be 
used for very flattened configurations, and it takes 100- 
200 iterations for the density and enthalpy distributions 
to converge. 
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Paper I mentions another numerical difficulty which 
has to do with the spike of the j(m ro ) curve near m ro = 1. 
The authors of Paper I have to truncate a small portion 
of the j(m ro ) curve in order to make the iteration con- 
verge. They also demonstrate that this truncation does 
not affect the inner structure of the star. It turns out 
that, for reasons still to be understood, the numerical 
instability associated with the curve only occurs 

for the most rapidly rotating models (i.e. those models 
where the pre-collapse white dwarfs have VL/Vt m = 1). 
Hence, the truncation is not necessary for all the other 
cases. 



0.25 




0.2 0.4 0.6 0.8 1 



FIG. 2. The values of /3 of the resulting neutron stars as a 
function of Q/Q m of the pre-collapse white dwarfs. 
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FIG. 3. Meridional density contours of the neutron star 
resulting from the AIC of a rigidly rotating O-Ne-Mg white 
dwarf with Q/Q m — 0.964. This neutron star has — 0.255. 
The contours in the upper graph denote, from inward to out- 
ward, p/p c = 10~ 4 , 10~ 6 , 10" 7 , 10" 8 , 10 -9 and 0. 
The contours in the lower graph denote, from inward to out- 
ward, p/p c = 0.8, 0.6, 0.4, 0.2, 0.1, 10" 2 , 10" 3 and 10" 4 . 
The central density of the star is p c = 3.3 x 10 14 g cm -3 . 



As in the case of pre-collapse white dwarfs, we measure 
the accuracy of the equilibrium models by the quantity e 
defined in Eq. (||). The models computed in this subsec- 
tion have e ranges from about 10~ 6 (for slowly rotating 
stars) to 10~ 4 (for rapidly rotating stars). We construct 
a number of neutron star models resulting from the col- 
lapse of the three sequences of pre-collapse white dwarfs 
in the previous subsection. The ratio of rotational ki- 
netic energy to gravitational potential energy, (3, of the 
neutron stars is plotted in Fig. || as a function of f2/Q m 
of the pre-collapse white dwarf. The values of /3 for all 
the neutron star models are smaller than 0.27, the critical 
value of P for the dynamical instability of rigidly rotating 
Maclaurin spheroids. 



The structure of the neutron stars with (3 > 0.1 are all 
similar: they contain a high-density central core of size 
about 20 km, surrounded by a low-density torus-like en- 
velope. The size of the envelope depends on the amount 
of rotation of the star, which can be measured by (3. The 
size ranges from 100 km (for (3 ~ 0.1) to over 500 km 
(for (3 Z 0.2). Figure || shows the density contours of a 
typical model. This model corresponds to the collapse 
of an O-Ne-Mg white dwarf with Q/Q m = 0.964. The 
resulting neutron star has (3 — 0.255. The envelope ex- 
tends to about 1530 km in this case. As a comparison, 
the equatorial radius of the pre-collapse white dwarf is 
2156 km. 
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FIG. 4. The distribution of rotational angular velocity Q 
as a function of zu for zu < 60 km. These are models for se- 
quence III with /3's of the resulting neutron stars equal to (a) 
0.0106, (b) 0.0555, (c) 0.0860, (d) 0.124, (e) 0.169, (f) 0.208, 
(g) 0.238 and (h) 0.261. The equatorial radii of the neutron 
stars in cases (a) and (b) are smaller than 60 km, and their 
frequency curves terminate at their equatorial radii. 




03 (km) 

FIG. 5. The cylindrical mass fraction function of 

zu for neutron star models in Fig. [|. The curves and their 
corresponding models are identified by the degree of central 
condensation: the higher the degree of central condensation, 
the lower the value of fi. 




m (km) 

FIG. 6. The value of fin, as a function of zu for the neutron 
star models in Fig. ^. The curves for models (a)-(g) terminate 
at the equatorial radii of the stars. 

Figure [| shows the rotational angular velocity Q as a 
function of radius w for neutron star models correspond- 
ing to the collapse of sequence III white dwarfs. The 
cases for sequences I and II are similar. We see that stars 
with small (3 show strong differential rotation. However, 
the rotation in the core region (m ^ 20 km) becomes 
more and more rigidly rotating as ft increases. The most 
rapidly rotating case (ft = 0.261) has £7 « 4400 rad cm -1 
in the core. This corresponds to a rotation period of 
1.4 ms, slightly less than the period of the fastest ob- 
served millisecond pulsar (1.56 ms). Further analysis 
reveals that the rotation curve in the envelope region 
roughly follows the Kepler law O oc tu~ 3 ^ 2 . 

Figure [5] shows the cylindrical mass fraction m CT as 
a function of tu for the same models as in Fig. ||. As 
expected, the degree of central concentration decreases 
with increasing ft. However, more than 80% of the mass 
is still concentrated inside a radius w — 30 km, even for 
the most rapidly rotating case. The collapsed object can 
be regarded as a rotating neutron star surrounded by an 
accretion torus. 

Numerous numerical studies demonstrate that the 
quantity ft is an important parameter for the dynami- 
cal stability of a rotating star. It is then useful to define 
a function as 

_ duu> tu' dz> [w'£l(LO>)Yp{w', z>) 
/ro | dw> m> dz' p{w>, js0$(tj7', z')\ ' ( > 

which is a measure of ft for the material inside a cylin- 
der of radius tu from the rotation axis. Figure || plots 
/3 ro as a function of tu for the same neutron star mod- 
els as in Fig. ||. We see that the curves level off when 
tu > 20 — 100 km for all rapidly rotating models, suggest- 
ing that the material outside 100 km is probably unim- 
portant for dynamical stability. It should also be no- 
ticed that the major contribution to ft is from the region 
10 km < tu < 50 km. Hence we expect that the material 
in this region plays an important role on the dynamical 
stability of the star. 
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III. STABILITY OF THE COLLAPSED OBJECTS 



dtC a + v b V b £ a - £ b V b v a = Sv a 



(18) 



In this Section, we study the dynamica l sta bilities of 



the neutron star models computed in Sec. II B using the 
technique of the linear stability analysis developed by 
Toman et al p9[ . This technique is briefly reviewed 
We then report the stability results in 



in Sec. 


[II A 


Sec. HI B 





The Eulerian change of the gravitational potential <5$ 
satisfies the Poisson equation 



V a V a (5$ = AttG8 P . 



(19) 



We find it useful to introduce a quantity Sh = SP/p, 
which is related to Sp by 



A. Linear stability analysis 

The motion of fluid inside the star is described by the 
hydrodynamical equations: 



d tP + V a {pv a ) = , 

, V a P 

d t v a + v b V b v a = V a $ 

P 

V a V a $ = AnGp , 



(9) 
(10) 

(11) 



where the summation convention is assumed, and V a 
denotes the covariant derivative compatible with three- 
dimensional flat-space metric. To study the stability, we 
perturb the density p and velocity v a away from their 
equilibrium values by small quantities: 



p(x b ,t) = p (x b )+6p(x b ,t) , 
v a (x b , t) = vjQ,{m)e% + Sv a (x b , t) , 



(12) 
(13) 



where is the unit vector along the azimuthal direction. 
The Lagrangian pressure perturbation AP is related to 
the Lagrangian density perturbation Ap by 



AP = lp -Ap , 
P 



(14) 



where for simplicity, the subscript "0" is suppressed, and 
hereafter in this Section, p and P denote the equilibrium 
density and pressure respectively. The quantity 



7 P 



dlogP 
dlogp 



(15) 



is the adiabatic index for pulsation. The relation between 
the Eulerian perturbations SP and bp can be easily de- 
duced from the transformation between the Lagrangian 
and Eulerian perturbations. The result is 



S p = j p -Sp+ ( 1 ) 

P V7cq 



where 



7eq 



dlogP 
dlogp 



(16) 



(17) 



cq 



is the adiabatic index computed from the equilibrium 
EOS. The Langrangian displacement £ a satisfies the 
equation 



6h = 7p ^5p+(^ -l)CVah. 

P \7cq 



(20) 



In the region where 7 P = 7 oq , Sh is the Eulerian change 
of the enthalpy. 

If the system is unstable, the perturbed quantities will 
grow in time. Instead of solving the fully non-linear equa- 
tions (pl)-([ll"|), however, Toman et al [29| develop a more 
efficient approach: expand Eqs. (||)-(11) to linear order 
of the perturbations and evolve the linearized equations. 

Consider the angular Fourier decomposition of any per- 
turbed quantity 5q: 



5q(x b ,t) = ^2 &qrn(™,z,t)t 



(21) 



It can be easily proved that each m-mode decouples in the 
linearized hydrodynamical equations because of the ax- 
isymmetry of the equilibrium configuration. In addition, 
the fact that the equilibrium configuration is symmet- 
ric under reflection about the equatorial plane (z — > —z) 
implies the modes with even and odd parity under the 
transformation z — > — z also decouple. Hence each m- 
mode with a definite parity can be evolved separately and 
the 3+1 simulation is reduced to a 2+1 simulation, which 
saves a lot of computation time. Hereafter, all perturbed 
quantities will be assumed to have angular dependence 

e irrup 

In Ref . |2£j , Toman et al choose to evolve the variables 
Sp and 8v a . However, we find it more convenient and 
numerically stable in our case to evolve the variables Sh 
and Sp a = pSv a . The reason being that the simulations 
are performed on a discrete grid, and it is preferable to 
use variables that change smoothly to ensure accuracy. 
However, the background density p decreases abruptly 
outside the core region, and the perturbation Sp is ex- 
pected to behave similarly. On the other hand, h, and 
presumably Sh, change much more smoothly even near 
the boundary of the star. In the case where 7 P + 1 7 oq , 
we also need to evolve the scalar function ?y = £ b Vbh. 
In terms of the new variables, the linearized equations 
become 

P. 



dtSh = —imflSh — 7 P — V a i5p a 

t (i-.)n.,, 

V7cq / P 

d t S P a = ~5p b V b v a - v b V b Sp a - P V a Sh - pV a <5$ 



(22) 
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lp_ 

7eq 



1 



{5h + r])W a h , 



d t r] 



Sp a 

-imfli] H V a /i 

P 



AnGp 2 
1 P P 



(23) 
(24) 
(25) 



It follows from Eqs. (gg)-(g5|) that if {5h,6p a ) is a solu- 
tion for an ro-mode, the complex conjugate (Sh* , Sp a *) is 
a solution for the -m-mode. We can then define the phys- 
ical "enthalpy" perturbation 5h = Sh + Sh* , and similarly 
for the physical density bp and velocity Sv a perturbations 
of an m-mode. 

We use a uniform cylindrical grid to perform the sim- 
ulations. We have checked that the code is able to repro- 
duce the results in Ref. [^9|. However, unlike the case in 
Rcf. [^9|, the collapsed objects studied here have a large 
envelope extending beyond 1000 km when the stars un- 
der consideration are rapidly rotating. This numerical 
difficulty can be circumvented by a suitable truncation 
scheme. 



As pointed out in Sec. [IB , we expect the outer enve- 
lope will not influence the dynamical stability in any sig- 
nificant way. Hence it is necessary to evolve the pertur- 
bations only in the dynamically interesting region. This 
is done by introducing a radius R m and a minimum den- 
sity p m in = p(Rm,Q). The perturbations are set to zero 
wherever the equilibrium density p(tu, z) < p m i n . Ifi? m is 
sufficiently large, increasing its value will not change the 
evolution result. We find that a value of R m » 200 km 
is needed to ensure that the results converge, and we use 
a cylindrical grid with 400 x 400 grid points to achieve a 
resolution of 0.5 km. 

In general, the two adiabatic indices j p and 7 oq are 
not equal. They coincide only if the pulsation timescale 
is much longer than all the reaction timescales for the 
different species of particles in the fluid to achieve equi- 
librium. This is the case for densities below neutron drip 
(p < 4 x 10 11 g cm -3 ) and above about 10 13 g cm" 3 . 
However, in the density range 4 x 10 11 g cm" 3 < p < 
10 13 g cm" 3 , the matter is a mixture of electrons, neu- 
trons and nuclei in equilibrium. Some of the reactions re- 
quired to achieve equilibrium involve weak interactions, 
which have timescales much longer than the pulsation 
timescale. Hence equilibrium is not achieved during pul- 
sation, and 7 P 7^ 7 cq in that density range Most 
people studying neutron star pulsations neglect the dif- 
ference of 7 P and 7 eq and use 7 oq in their calculations. It 
has been demonstrated (see e.g. [^3| ) that this treatment 
has no significant effect on the final result, because the 
matter in that density range occupies only a tiny fraction 
of neutron star. However, it may have an important ef- 
fect on the stability of the new-born neutron stars studied 
here. The reason is that th e dy namically important re- 
gion, as pointed out in Sec. |IB|, is 10 km < w < 50 km. 
This region contains a significant amount of matter in 
that density range (see Fig. ||). Our numerical simula- 



tions indicate that this is indeed the case. The critical 
value Pd for the dynamical instability drops from about 
0.25 to 0.23 if 7 oq is used for the adiabatic index of pul- 
sation. 




5 7 9 11 13 15 

log 10 [p (g cm" 3 )] 

FIG. 7. The values of 7 cq (solid line) and 7 P (dashed line) 
as a function of log 10 p. The two curves coincide when p < pd 
and p > p e . 

The appropriate 7 P remains roughly constant from 
the density of neutron drip pd ' 1 " 1 ' 
to the density above which 7 oq 

g cm -1 " 
7 p which mimics the curve in Refs. [ |30| , [31| and which is 
compatible with the EOS used here, we take 7 P in the 
density range pd < p < p e to be (also see Fig. ^) 



x 10^ g 
around p e 



10 13 g cm -3 ||30|,pli. To provide a reasonable value of 



lp{p) = leq(Pd) + [leq(Pe) ~ Jcq(pd)} 



log 2 (p/pd) 
^g 2 {p e /pd) 



(26) 



Under some circumstances it is possible to have a re- 
gion of the star where the mode is stationary in the 
fluid's co-rotating frame. In this case, we should use 7 eq 
for the adiabatic index of pulsation in the region where 
\uj'\ = \uj+mfl\ <C 2n/t r . Here u> is the angular frequency 
of an m-mode that has dependence exp[i(ujt-\-mip)] in the 
inertial frame; to' = lo + mQ is the angular frequency of 
the mode in the fluid's co-rotating frame; and t r w 1 s is 
the timescale for different species of particles to achieve 
/3-equilibrium in the density range pd < p < p e - It 
turns out (see the next subsection) that rapidly rotat- 
ing neutron stars have an unstable bar mode (m = 2) 
with oj « —3000 rad s _1 . There is indeed a radius at 
which to' = 0. This radius is at w — vj c w 40 km for 
stars with (3 > 0.23. The density on the equator of the 
stars is p « 10 12 g cm" 3 , well within the questionable 
density range. However, the width of this "co-rotating 
region" which satisfies \oj'\ < 2ir/t r is 



Aw 



2-ir/t r 



|9„n(w c )| 



0.1 km 



The material in the region contains only 10 4 of total 
mass and angular momentum of the star. Hence this 
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thin co-rotating layer is not expected to have a signifi- 
cant influence on the overall stability of the stars. 



B. Results 

We perform a number of simulations on neutron star 
models computed by the method described in Sec. II B. 
The simulations are terminated either when an instabil- 
ity is fully developed or when the simulation time reaches 
60 ms, corresponding to 40 rotation periods of the ma- 
terial at the center of the star. We regard a star as dy- 
namically unstable if the density perturbation shows an 
evidence of exponential growth and increases its ampli- 
tude by at least a factor of fifteen by the end of the sim- 
ulation. In our simulations, no instability is observed for 
neutron star models in sequences I and II. A bar-mode 
(rn = 2) instability develops for sequence III models when 
the star's (3 is greater than a critical value (3d ~ 0.25. The 
unstable mode has even parity under reflection about the 
equatorial plane. This (3d is slightly less than the criti- 
cal value 0.27 for the Maclaurin spheroids. It should be 
pointed out that all the stars in sequences I and II have 
(3's smaller than this fid- Hence we believe that they are 
stable simply because their /3's are not high enough. 

Some other simulations fi|ll| show that in the cases 
where (3d < 0.27, the instability is dominated by the 
m = 1 mode for stars with f3 close to (3d- However, we 
do not observe any sign of an unstable m = 1 mode in 
our case. We also performed simulations using 7 eq (the 
solid curve in Fig. |7]) as the adiabatic index for pulsation 
instead of j p (the dashed curve in Fig. 0). We find that 
(3d drops to about 0.23, showing that matter in the den- 
sity region pd < p < p e plays an important role on the 
instability. 



1000 - 




FIG. 8. The relative amplitude a r as a function of time. 
The equilibrium star has (3 = 0.261, the most rapidly rotat- 
ing model. 



To visualize the instability, we define an amplitude 

■j\s P \ 2 d 3 x^ 1/2 



a 



J p 2 d 3 a 



(27) 



for the density perturbation. Since we evolve the pertur- 
bations using linearized equations, it is more convenient 
to work with the relative amplitude a r : 



a r (t) = a(t)/a(0) . 



(28) 



This relative amplitude is defined so that it is equal to 
one at t = 0. Figure || shows the time evolution of a r 
for the most rapidly rotating star (J3 = 0.261). We see 
that after about 10 ms, an instability develops and a r 
grows exponentially. The e-folding time of the growth r 
is found, by least square fit, to be 7.8 ms. 
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FIG. 9. The time evolution of the angular frequency 0^22 
for the most rapidly rotating star {fi — 0.261). 

The unstable mode can also be characterized by a com- 
plex angular frequency defined as 



D 



W22 



22 



1D0 



where 



D 22 = / r 2 5pY* 2 d 3 x 



(29) 



(30) 



is the mass quadrupole moment, and the spherical har- 
monic function 

Y 22 = l.L^si^ee 2 ** . 
4 V 2tt 

The time derivative of D22 is evaluated by the for- 
mula l3l 



D22 = I Sp a V a {r 2 Y* 2 )d 3 x-2i I Qr 2 6pY* 2 d 3 x , (31) 

where we have used the continuity equation (^) and in- 
tegrated by parts. 

Let w be the complex frequency of the most unsta- 
ble mode. The e-folding time is related to the imagi- 
nary part of u> by Im(aj) = — 1/t. At late time, the 
density perturbation is dominated by the most unstable 
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mode, which means that both 8p and D22 go approxi- 
mately as exp(iwi). Hence 0J22 ~ w. Fig. || plots LO22 as 
a function of time for the evolution of the most rapidly 
rotating star. We see that at late time, W22 is approx- 
imately a constant, indicating that the perturbation is 
indeed dominated by the most unstable mode. The fre- 
quency of the unstable mode is then determined to be 
u> w (—2890 — 130i) rad s _1 . Note that the imaginary 
part agrees with the e-folding time determined above. 




10 20 30 40 50 60 

JjS (km) 

FIG. 10. The magnitude of the density perturbation \5p\ 
of the unstable bar mode of the most rapidly rotating star 
(/3 — 0.261) on the equatorial plane. The magnitude is nor- 
malized so that the maximum value is one. 

Figure [lO] shows the magnitude of the density pertur- 
bation \Sp\ of the unstable bar mode of the most rapidly 
rotating star {(3 = 0.261) on the equatorial plane. We 
see that \Sp\ has a peak at w ~ 20 km, which is in the 
transition region between the neutron star core and the 
tenuous outer layers (see Fig. ||). There is a small, sec- 
ondary peak at vo = 44 km, which is the corotation radius 
at which Re(u>) + 2f2 = for this neutron star. This sec- 
ondary peak is caused by the resonant response of the 
fluid being driven by the mode corotating with it (see 
Appendix H). 



10 




o o o o o o o 
UJ ■+ oj cn ^ to 

I I I 

X (km) 

FIG. 11. The eigenfunction of the physical density pertur- 
bation Sp on the equatorial plane for the bar mode of the 
most rapidly rotating star (/3 = 0.261). The density pertur- 
bations are normalized so that the maximum value is one. 
Only positive density regions of the eigenfunction are shown. 
The negative structure of the eigenfunction can be inferred 
from the sinusoidal structure of the eigenfunction. The con- 
tour levels are, from inward to outward, 0.8, 0.6, 0.4, 0.2, 0.1 
and 0.01. The small arcs inside the 0.1 contours are addi- 
tional contours of 0.2, corresponding to the secondary peak 
in Fig. [o| 
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FIG. 12. The eigenfunction of the physical velocity pertur- 
bation Sv a on the equatorial plane for the bar mode of the 
most rapidly rotating star (f3 = 0.261). 



Figures |ll| and |12| show the eigenfunctions of the phys- 
ical perturbations Sp and Sv a on the equatorial plane. 
Note that our grid extends out to 200 km from the 
center, but the dynamically interesting region is con- 
centrated within 60 km from the center. Since the 
time dependence of the perturbations go as exp[i(u>t + 
rrup)\, Re(w) < means the pattern rotates in prograde 
(counter-clockwise) direction. The density perturbation 
is bar-like in the inner region and becomes trailing spi- 
rals in the outer region. Similar structure is also observed 
in other numerical simulations on the bar-mode instabil- 
ity ^QIQUPQ. The secondary peak of Sp appears 
as two small arcs in Fig. [Il] just inside the 0.1 contours. 
Figure |l2| shows that Sv a is almost parallel to the (p di- 
rection at the corotation radius, which is also a result of 
resonance (see Appendix ^) . Since Sv a changes abruptly 
near the corotation radius, it is very possible that shocks 
will develop there when the perturbations become large. 
This might have significant influence on the non-linear 
evolution of the bar mode. 

TABLE II. The oscillation frequency / and e-folding time 
r of the most unstable bar mode for several unstable neutron 
stars. Here Q is the rotational frequency of the pre-collapse 
white dwarf, and £l m is the maximum frequency of the white 
dwarf in Sequence III. 



(n/n m )w£) 





Re(-w) 
rad s 


/ 

Hz 


T 

ms 


0.934 


0.251 


2800 


445 


20 


0.964 


0.255 


2850 


450 


12 


0.989 


0.258 


2850 


450 


8.9 


1.000 


0.261 


2890 


460 


7.8 



11 



The eigenfunctions of the most unstable bar mode for 
the other unstable equilibrium neutron stars are similar 
to those displayed above. Table |l| summarizes the oscil- 
lation frequencies [/ = |Re(u;)|/(27r)] and e-folding time 
r of the unstable models we have studied. The table also 
shows the ratio of the rotational frequency of the pre- 
collapse white dwarfs to the maximum frequency fl m of 
the white dwarf in the sequence. We find that the os- 
cillation frequencies are almost the same (« 450 Hz) for 
all the cases. We do not observe any instability in our 
simulations for stars with (3 < 0.241. Hence we conclude 
that (3d is somewhere between 0.241 and 0.251, and the 
pre-collapse white dwarf has to have Q > 0.930m in order 
for the collapsed star to develop a dynamical instability. 



IV. GRAVITATIONAL RADIATION 

In this section, we estimate the strength of the grav- 
itational radiation emitted by neutron stars undergoing 
a dynamical instability. We also estimate the signal to 
noise ratio and discuss the detectability of these sources. 

The rms amplitude of a gravitational wave strain, h(t), 
depends on the orientation of the source and its location 
on the detector's sky. When averaged over these angles, 
its value is given by pi] 



h*(t) = ~(hl(t) + hl(t)) 



(32) 



where h + (t) and h x (t) are the rms amplitudes of the plus 
and cross polarizations of the wave respectively, and (...) 
denotes an average over the orientation of the source and 
its location on the detector's sky. 

In the presence of perturbations, the density and ve- 
locity of fluid inside the star become 

oo 

p{x,t) = p (w, z) + ^2 Sp m (x,t) (33) 

m— — oo 

oo 

v(x, t) = wQ,(vo)e,f, + y~] Sv m (x,t) , (34) 



m— — oo 



where the perturbation functions 5p m and Sv m have an- 
gular dependence e" n(p . The amplitude of the gravita- 
tional waves produced by time varying mass and current 
multipole moments can be derived from Ref. |f37f . The 
result is 



1 



l 



1=2 m=- 



AG 2 
hd 21 



\D 



CO |2 
lm I 



|sfl 2 



(35) 



where D is the distance between the source and detector; 
c is the speed of light, and 



N, = 



4tt(Z + 1)(Z + 2) 
Z(/-l)[(2l + l)!!] a 



(36) 



D 



(i) 



S, 



a 1 



dtl D lm ; 



(37) 
(38) 



lm dt l ' 

For a Newtonian source, the mass moments Di m and cur- 
rent moments Si m are given by 

/),,„: / pr')Y„, d 3 X 



(39) 



S, 



I 



l + l 



r l pv 



1 lm 



d 3 x , 



(40) 



where 



x x VYi m / + 1) are the magnetic 
type vector spherical harmonics. The functions Di m and 
Sim have the property that D* m = (— l) m Z?/_ m and 
Hence it is sufficient to consider 



lm 



(-l) m Si. 



only positive values of m and Eq. (135]) becomes 



1 



8G 2 



D 2 5c 2/+4 
1=2 m=0 



N, 



\D 



(0|2 
lm\ 



IS' 



(0|2 

lm I 



(41) 



The energy and angular momentum carried by the 
gravitational waves can also be derived from |57). The 
result is 



^=EE4^ 



1=2 m=0 
oo / 



\D 



(i+i) 



\S, 



a+i) 



G 



1=2 m=0 



(I)* (l+l) q (l)* q (l+l) 

^Im ^Im "T" °lm °lm 



(42) 
(43) 



where the overline denotes time average over several pe- 
riods. 

When a neutron star develops a dynamical instabil- 
ity and the bar mode (m = 2) is the only unstable 
mode, the values of h, E and J will be dominated by 
the term involving D22 ■ Since the unstable bar mode has 
even parity under reflection about the equatorial plane, 
D32 = S22 = and the next leading term will involve S32 
and Z?42. These terms are expected to be smaller than 
the D22 term by a factor of (v/c) 4 for E and J, and a 
factor of {v/c) 2 for h. In our models, v/c < 0.28, so the 
contribution of higher order mass and current multipole 
moments are small and will be neglected. 

Strictly speaking, the above analysis only applies when 
the amplitudes of the perturbations are small. When the 
amplitudes are large, however, the fluid motion does not 
separate neatly into decoupled Fourier components, so 
all Dim and Si m will contribute. However, it is expected 
that the D22 term will still be the most important term. 
Since the detailed non-linear evolution of the dynami- 
cal instability is not known, the aim of this section is to 
provide an order of magnitude estimate of the gravita- 
tional radiation from these sources. Hence we shall only 
consider the effect of the mass quadrupole moment and 
assume D22 can be approximated by the bar-mode eigen- 
functions computed in Sec. Ill B . In this approximation, 
Eqs. ©-© become 
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E = 

j = 



5c 4 D 
1024tt 7 G 

75c 5 
1024tt 6 G 
75c 5 



7T 

15 ; 



/ 6 ^22| 2 



22 



(44) 
(45) 
(46) 



where / = |Re(w)|/(27r) is the oscillation frequency of 
the bar mode. 

Substituting the bar-mode eigenfunctions (from 
Sec. |ffl B[) into Eq. @, we find that 



\D 



22 



a 8 x 10 45 g cm 



(47) 



for all the unstable models we have studied. Here a is the 
amplitude of the bar mode defined in Eq. (p^j) . The mass 
quadrupole moment D22 has a time dependence exp(icot), 
where 10 is the angular frequency of the mode. Hence the 
time derivative -D 2 2 = {i<-^) 1 D22 and we obtain 



_ , ~ 3 ,'20 Mpc\ 
ft«a7x 10 2J 1 ' 



D 



E^a 2 9 x 10 52 erg s _1 : 
J«a 2 6 x 10 49 g cm 2 s" 



5 



(48) 

(49) 
(50) 



The signal to noise ratio of these sources depends on 
the detailed evolution of the bar mode when the density 
perturbation reaches a large amplitude and non-linear ef- 
fects take over. Recently, New, Centrella and Tohlinc |Tl| ] 
and Brown ]3qj perform long-duration simulations of the 
bar-mode instability. They find that the mode saturates 
when the density perturbation is comparable to the equi- 
librium density, and the mode pattern persists, giving a 
long-lived gravitational wave signal. Here we assume that 
this is the case, and that the mode dies out only after 
a substantial amount of angular momentum is removed 
from the system by gravitational radiation. We then fol- 
low the method described in Refs. |38],[39| to estimate the 
signal to noise ratio. 

In the stationary phrase approximation, the gravita- 
tional wave in the frequency domain h(f) is related to 
h(t) by 



h 2 (t) = \h(f)\ 2 
Combining Eqs. ([I3J) 

\Hf)\ 2 



'1 
dt 

|) and (fl|), we obtain 
G j 



c A 57r/|/|£> 2 
The signal to noise ratio is given by 



^ =2 



\kf)\ 2 

Sh(f) 



df , 



(51) 



(52) 



(53) 



where Sh(f) is the spectral density of the detector's noise. 
If we assume that the oscillation frequency remains con- 
stant in the entire evolution, we obtain |4Cf| 



S_ 

N 



2G AJ 



DV 5nc?fS h (f) ' 



(54) 



where AJ is the total amount of angular momentum 
emitted by gravitational waves. To estimate AJ, we as- 
sume that the mode dies out when the angular momen- 
tum of the star decreases to Jj « 3.3 x 10 49 g cm 2 s -1 , 
which is the angular momentum of the marginally bar- 
unstable star. Then we have AJ < 5 x 10 48 g cm 2 s _1 
for all the unstable stars, and the signal to noise ratio for 
LIGO-II broad-band interferometers |4lJ is 



S_ 

N 



■— = 15 



20 Mpc\ 



D 



AJ 



5 x 10 48 cgs 



1/2 



( f 



V450 Hz 



-1/2 



VMf) 



2 x 10~ 24 Hz 



-1/2 



(55) 



The timescale of the gravitational wave emission can 
be estimated by the equation 



tgw 



AJ 
j 

- / a s 
78 Ut 



AJ 



5 x 10 48 cgs 



(56) 



where a s is the amplitude a of the density perturbation 
at which the mode saturates. We have used Eqs. ( |46| ) 
and (|47]) to calculate the numerical value in the last equa- 
tion. 

The detectability of this type of sources also depends 
on the event rate. The event rate for the AIC in a 
galaxy is estimated to be between 10 -5 and 10 -8 per 
year [^2|,Q. Of all the AIC events, only those corre- 
sponding to the collapse of rapidly rotating O-Ne-Mg 
white dwarfs can end up in the bar-mode instability, and 
the fraction of which is unknown. If a signal to noise 
ratio of 5 is required to detect the source, an event rate 
of at least 10 _6 /galaxy/year is required for such a source 
to occur at a detectable distance per year. Hence these 
sources will not be promising for LIGO II if the event 
rate is much less than 10 -6 per year per galaxy. 

The event rate of the core collapse of massive stars is 
much higher than that of the AIC. The structure of a pre- 
supernova core is very similar to that of a pre-collapse 
white dwarf, so our results might be applicable to the 
neutron stars produced by the core collapse. If the core 
is rapidly rotating, the resulting neutron star might be 
able to develop a bar-mode instability. If a significant 
fraction of the pre-supernova cores are rapidly rotating, 
the chance of detecting the gravitational radiation from 
the bar-mode instability might be much higher than ex- 
pected. 
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V. MAGNETIC FIELD EFFECTS 



As mentioned in Sec. II B, a new-born hot proto- 



neutron star is dynamically stable because its f3 is too 
small. It takes about 20 s for the proto-neutron star to 
cool down and evolve into a cold neutron star, which may 
have high enough (5 to trigger a dynamical instability. 
The proto-neutro n star s, as well as the cold neutron stars 
computed in Sec. II B, show strong differential rotation 
(Paper I) . This differential rotation will cause a frozen-in 
magnetic field to wind up, creating strong toroidal fields. 
This process will result in a re-distribution of angular 
momentum and destroy the differential rotation. If the 
timescale of this magnetic braking is shorter than the 
cooling timescale, the star may not be able to develop 
the dynamical instability discussed in Sees. [II and [TV]. 
In this Section, we estimate the timescale of this mag- 
netic braking. 

In the ideal magnetohydrodynamics limit, the mag- 
netic field lines are frozen into the moving fluid. The 
evolution of magnetic field B is governed by the induc- 
tion equation 



dB a 
~~dt~ 



v b V b B a - B b V b v a = -B a V b v b 



(57) 



In our equilibrium models, v b = vjVl(vj)e b i . Hence 
V b v b = and Eq. (|57|) becomes 



dB a 
dt 



B b V b v a , 



(58) 



where d/dt — d/dt + v b W b is the time derivative in the 
fluid's co-moving frame. Eq. (p8 [ ) c an be integrated ana- 
lytically (see e.g. Appendix B of fl44|). The magnetic field 
B^(x,t) at the position a; of a fluid element at time t is 
related to the magnetic field B k (x ,t ) at the position 
x Q of the same fluid at time t by 



B j {x,t) 



(59) 



where dx^ jdx\ is the coordinate strain between to and t. 
With v a = wH,(zu)e'i, it is easy to show that the induced 
magnetic field has components only in the direction. 
Its magnitude Bi, after a time t, is easy to compute from 
Eq. (M). The result is 



Bi(t) = B tvu\d^n\ 



(60) 



where Bo is the component of magnetic field in the e ro 
direction. The induced magnetic field will significantly 
change the equilibrium velocity field when the energy 
density of magnetic field €b — Bf/(8n) is comparable 
to the rotational kinetic energy density €r = pw 2 ft 2 /2. 
This will occur in a timescale tb set by es = £r- Using 
Eq. (pOj), we obtain 



TB 



M Bo 



L 



(61) 



where L = f2/|d ra f2| is the length scale of differential ro- 
tation, and va = B / ' ^/Airp is the speed of Alfven waves. 

Observational data suggest that the magnetic fields of 
most white dwarfs are smaller than 10 5 G, although a 
small fraction of "magnetic white dwarfs" can have fields 
in the range 10 6 — 10 9 G. Assuming flux conservation, 
the magnetic fields of the hot proto-neutron stars just af- 
ter collapse would be Bo ~ 10 9 G for those 10 5 G white 
dwarfs. Using the angular velocity distribution in Pa- 
per I for the hot proto-neutron star, we find that the 
magnetic timescale in the dynamically important region 
(w < 100 km) is 



t b « 10 4 s 



10 9 G 
Bo 



(62) 



which is much longer than the neutrino cooling timescale 
(~ 20 s) . Hence the angular momentum transport caused 
by the magnetic field is negligible during the cooling pe- 
riod. The magnetic timescale for the cold neutron stars 
can be calculated from the angular frequency distribution 
computed in Sec. II B. We find that tb for the cold mod- 

which is still 



els is about half of that given by Eq. fl62|) 
much longer than the timescale of gravitational waves 
tqw calculated in the previous Section. The instabil- 
ity results presented in the previous two Sections remain 
unchanged unless the neutron star's initial magnetic field 
Bq is greater than 10 12 G. In that detailed mag- 

netohydrodynamical simulation has to be carried out to 
compute the angular momentum transport. 

The magnetic timescale for these nascent neutron stars 
is significantly different from that estimated by Baum- 
garte, Shapiro and Shibata |4|] and Shapiro jig ]. They 
consider differentially rotating "hypermassive" neutron 
stars, which could be the remnants of the coalescence of 
binary neutron stars. Those neutron stars are very mas- 
sive (M ~ 3M Q ) and have much higher densities than 
the new-born neutron stars studied in this paper. They 
also use a seed magnetic field of strength Bo ~ 10 12 G, 
which is much larger than our estimate. These two dif- 
ferences combined make our magnetic braking timescale 
two orders of magnitude larger than theirs. It should be 
noted that it is the magnetic field just after the collapse 
that is relevant to our analysis here. The strong differen- 
tial rotation of the neutron star will eventually generate 
a very strong toroidal field (Bi ~ 10 16 G) and destroy 
the differential rotation. The final state of the neutron 
star will be in rigid rotation, and its magnetic field will 
be completely different from the initial field. For this rea- 
son, the field strength B ~ 10 12 G observed in a typical 
pulsar is probably not relevant here. 



VI. SUMMARY AND DISCUSSION 



We have applied linear stability analysis to study the 
dynamical stability of new-born neutron stars formed by 
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AIC. We find that a neutron star has a dynamically un- 
stable bar mode if its (3 is greater than the critical value 
(3d ~ 0.25. In order for the neutron star to have j3 > (3d, 
the pre-collapse white dwarf must be composed of oxy- 
gen, neon, magnesium and have a rotational angular fre- 
quency f2 > 3.4 rad s _1 , corresponding to 93% of the 
maximum rotational frequency the white dwarf can have 
without mass shedding. 

The eigenfunction of the most unstable bar mode is 
concentrated within a radius w < 60 km. The oscillation 
frequency of the mode is / ~ 450 Hz. When the ampli- 
tude of the mode is small, it grows exponentially with 
an e-folding time t w 8 ms for the most rapidly rotating 
star ((3 = 0.261), which is about 5.5 rotation periods at 
the center of the star. 

The signal to noise ratio of the gravitational waves 
emitted by this instability is estimated to be 15 for 
LIGO-II broad-band interferometers if the source is lo- 
cated in the Virgo cluster of galaxies (D = 20 Mpc). 
The detectability of these sources also depends on the 
event rate. The event rate of AIC is between 10~ 5 and 
10 _8 /galaxy/year. Only those AIC events corresponding 
to the collapse of rapidly rotating O-Ne-Mg white dwarfs 
can end up in the bar-mode instability. While it is likely 
that the white dwarfs would be spun up to rapidly ro- 
tation by the accretion gas prior to collapse JL7J, it is 
not clear how many of the AIC events are related to the 
O-Ne-Mg white dwarfs. If the event rate is less than 
10~ 6 /galaxy/year, it is not likely that LIGO II will de- 
tect these sources. However, the event rate of the core 
collapse of massive stars is much higher than that of the 
AIC. A bar- mode instability could develop for neutron 
stars formed from the collapse of rapidly rotating pre- 
supernova cores. If a significant fraction of the cores 
are rapidly rotating, the chance of detecting the grav- 
itational radiation from bar-mode instability would be 
much higher. 

If the pre-collapse white dwarf is differentially rotat- 
ing, the resulting neutron star can have a higher value 
of (3. The bar-mode instability is then expected to 
last for a longer time. However, any differential rota- 
tion will be destroyed by magnetic fields in a timescale 
tb ~ R/va, where R is the size of the white dwarf and 
va = B/^/Anp ~ By/R 3 /(3A1). For a massive white 
dwarf with M = 1.4M Q , 

/10 5 G\ / R \~ 1/2 
TB ~ 2 ^{—){lM^) < (63) 

which is much shorter than the accretion timescale. 
Hence rigid rotation is a good approximation for pre- 
collapse white dwarfs. 

The magnetic field of a neutron star is much stronger 
than that of a white dwarf. The timescale for a magnetic 
field to suppress differential rotation depends on the ini- 
tial magnetic field B of the proto-neutron star. If the 
magnetic field of the pre-collapse white dwarf is of or- 
der 10 5 G, the initial field will be B ~ 10 9 G according 



to conservation of magnetic flux. In this case, the mag- 
netic timescale is tb ~ 10 4 s. This timescale is much 
longer than the time required for a hot proto-neutron 
star to cool down and turn into a cold neutron star, 
and go through the whole dynamical instability phase. 
If So ~ 10 12 G, a significant amount of angular momen- 
tum transport will take place during the cooling phase. 
A detailed magnetohydrodynamical simulation has to be 
carried out to study the transport process in this case. 
However, such a strong initial magnetic field is possible 
only if the pre-collapse white dwarf has a magnetic field 
B > 10 s G. 

Finally, we want to point out that the collapse of white 
dwarfs will certainly produce asymmetric shocks and may 
eject a small portion of the mass. We expect that our 
neutron star models describe fairly well the inner cores 
of the stars but not the tenuous outer layers. Our stabil- 
ity results are sensitive to the region with tu 5; 100 km. 
The results could change considerably if the structure in 
this region is very different from that of our models. This 
issue will hopefully be resolved by the future full 3D AIC 
simulations. 
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APPENDIX A: RESONANCE AT THE 
COROTATION RADIUS 

We see from Figs. [l(]-|l2] that the bar-mode eigen- 
function has peciliar structures at the corotation radius 
(m sa 40 km) at which ui + 2f2 w 0. The density pertur- 
bation has a small peak and the velocity perturbation is 
almost parallel to the if direction. In this Appendix, we 
shall show that these are caused by the resonance of the 
fluid driven by the mode. 

For simplicity, we only consider the fluid's motion on 
the equatorial plane. Assume that the perturbations are 
dominated by a mode that goes as exp(zwi + irrup). We 
also assume that this mode is even under the reflection 
z — > — z. Hence we have £ z = and 8v z = 0. In cylindri- 
cal coordinates, the linearized Euler equation takes the 
form 

B r)P BP 
i(w + mtySv™ = 2Q6v v - — — + -^-5p - 3^5$ , 

P P 

%{w + mn)Sv v = -(wd m fl + 2n)Sv™ — + 5<$>) 

m \ p ) 

The density perturbation Sp is related to the pressure 
perturbation SP by 
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5p = ^-5P - ( 1 - ^ ) r d^p . 



> 



(A3) 



The tu-component of the Lagrangian displacement is 
given by 



5v° 



i(lu + mQ) 



(A4) 



Our numerical simulations show that SP is well- 
behaved and smooth near the corotation radius at which 
lo + mil 0. The perturbed gravitational potential 5$ 
is expected (and is confirmed by our numerical simula- 
tions) to be smooth since it depends on the overall dis- 
tribution o f th e density perturbation. We can then use 
Eqs. (|Al[)-(|A"4|) to express all the other perturbed quan- 
tities in terms of SP and Near the corotation radius, 
the expressions are: 



Sv^ - 

K 2 - 

B - 



w(k 2 + B) 



SP - 1 _ _!S 



7cq \ Sv^d^p 



SP 
P 

2il 



— + <5$ 



rud^fl 2 + 4ft 2 
dzaPdvjp ( \ Toq 
7 P 



im ( 5P 
— — + 5$ 
zu \ p 



1 



(A5) 

(A6) 

(A7) 
(A8) 
(A9) 



It follows from Eqs. (|A4|) and (|AJ) that if \5v™\ is not 
of order (oj + mVt) near the corotation radius, both |£ ro | 
and \Sp\ will be large. The large magnitude of the La- 
grangian displacement is caused by the fluid being driven 
in resonance by the mode. The large displacement of 
the fluid causes \Sp\ to be large due to the second term 
of Eq. (|A3|). This term arises becuase of the different 
compressibilities of stationary and oscillating fluid (i.e. 
7oq 7^ 7p)- In the case of the bar mode (m — 2), the 
corotation radius is located at vj c ss 40 km. The equi- 
librium density on the equator p(w c , 0) « 10 12 g cm~ 3 
and the stationary fluid is very compressible (7 oq ~ 0.7). 
The high compressibility of the stationary fluid make the 
background equilibrium density p drop rapidly as w in- 
creases, i.e. \d^p\ is large. The oscillating fluid is far less 
compressible (7 P = 1.35). As a result, when the oscillat- 
ing fluid moves to a new location, it does not expand or 
compress to an extent that can compensate for the dif- 
ference between the background densities at the old and 
new locations. Since both |£ ro | and |<9 ro p| are large, Sp 
is dominated by the second term of Eq. (A3) near the 
corotation radius. This explains the narrow s eco ndary 
peak of 8p seen in Fig. [h]. We see from Eq. (A7) that 
l^w^l ^ \5v™\ and Sv^ changes rapidly near the coro- 
tation radius, which explains the flow pattern seen in 

Fig. im 



[i 

[2 
[3; 

[5 
[6 

[-: 

[8 
[9 
[10; 
[11 
[12 
[13 
[14; 
[15; 
[16 
[17 

[is; 

[19 

[20 

[21 
[22 

[23 
[24 
[25' 

[26 

[27 

[28 

[29 



K. S. Thorne, in Black Holes and Relativistic Stars, ed. 
R. M. Wald (The University of Chicago Press, 1998). 
S. Chandrasekhar, Elliposidal Figures of Equilibrium 
(New Haven: Yale University Press, 1969). 
J. E. Tohline, R. H. Durisen and M. McCollough, Astro- 
phys. J., 298, 220, 1985. 

R. H. Durisen, R. A. Gingold, J. E. Tohline and A. P. 
Boss, Astrophys. J., 305, 281, 1986. 
H. A. Williams and J. E. Tohline, Astrophys. J., 334, 
449, 1988. 

J. L. Houser, J. M. Centrella and S. C. Smith, Phys. Rev. 
Lett., 72, 1314, 1994. 

S. Smith, J. L. Houser and J. M. Centrella, Astrophys. 
J., 458, 236 (1996). 

J. L. Houser and J. M. Centrella, Phys. Rev. D., 54, 7278 
(1996). 

B. K. Pickett, R. H. Durisen and R. H. Davis, Astrophys. 
J., 458, 714 (1996). 

J. L. Houser, Mon. Not. Roy. Astro. Soc, 209, 1069 
(1998). 

K. C. B. New, J. M. Centrella and J. E. Tohline, Phys. 
Rev. D, 62, 064019 (2000). 

J. N. Imamura and J. Toman, Astrophys. J, 444, 363 
(1995). 

J. M. Centrella, K. C. B. New, L. L. Lowe and J. D. 
Brown, Astrophys. J. Lett., 550, 193 (2001). 
N. Stergioulas and J. L. Friedman, Astrophys. J., 492, 
301 (1998). 

M. Shibata, T. W. Baumgarte and S. L. Shapiro, Astro- 
phys. J., 542, 453 (2000). 

M. Saijo, M. Shibata, T. W. Baumgarte and S. L. 
Shapiro, Astrophys. J., 548, 919 (2000). 

C. L. Fryer, D . E. Holz and S. A . Hughes, to appear in 
Astrophys. J., ^stro-ph/0106113| (2001). 



See e.g. R. Narayan and R. Popham, Astrophys. J. Lett., 
346, 25 (1989). 

Y. T. Liu and L. Lindblom, Mon. Not. Roy. Astro. Soc, 
324, 1063 (2001) (Paper I). 

K. Nomoto, in Proc. 13th Texas Symp. on Relativistic 
Astrophysics, edited by M. Ulmer, (World Scientific, Sin- 
gapore, 1987). 

K. Nomoto and Y. Kondo, Astrophys. J. 367, L19 (1991). 
E. Bravo and D. Garcfa-Senz, Mon. Not. Roy. Astro. Soc, 
307, 984 (1999). 

E. E. Salpeter, Astrophys. J., 134, 669 (1961). 

I. Hachisu, Astrophys. J. Lett. 61, 479 (1986). 

J.-L. Tassoul, Theory of Rotating Stars, Princeton Univ. 

Press (1978). 

H. A. Bethe and M. B. Johnson, Nucl. Phys. A, 230, 1 
(1974). 

G. Baym, H. A. Bethe and C. J. Pethick, Nucl. Phys. A, 
175, 225 (1971). 

S. Smith and J. M. Centrella, in Approaches to Numerical 
Relativity, edited by R. d'Inverno, Cambridge University 
Press, New York (1992). 

J. T. Toman, J. N. Imamura, B. K. Pickett and R. H. 
Durisen, Astrophys. J., 497, 370 (1998). 



1G 



[30] D. W. Meltzer and K. S. Thorne, Astrophys. J., 145, 514 [39] 
(1966). [40] 

[31] M. Colpi, S. L. Shapiro and S. A. Teukolsky, Astrophys. 

J., 339, 318 (1989). [41] 

[32] L. S. Finn and C. R. Evans, Astrophys. J., 588 (1990). 

[33] L. Lindblom and S. L. Detweiler, Astrophys. J., 53, 73 [42] 
(1983). 

[34] J. N. Imamura, R. H. Durisen and B. K. Pickett, Astro- 
phys. J., 528, 946 (2000). "[437 

[35] J. D. Brown, Phys. Rev. D., 62, 0004002 (2000). 

[36] See e.g. K. S. Thorne in Three Hundred Years of Gravita- [44] 
tion, edited by S. W. Hawking and W. Israel (Cambridge 
University Press 1987) . [45] 

[37] K. S. Thorne, Rev. Mod. Phys., 52, 299 (1980). 

[38] B. J. Owen, L. Lindblom, C. Cutler, B. F. Schutz, A. Vec- [46] 
chio and N. Andersson, Phys. Rev. D 58, 084020 (1998). 



B. J. Owen and L. Lindblom, |gr-qc/0111024| (2001). 
This formula was first derived by R. D. Blandford (un- 
published) . 

E. Gustafson, D. Shoemarker, K. Strain and R. Weiss, 
LIGO Document T990080-00-D. 

V. Kalogera, to appear in the proceedings of the 



3rd Amaldi Conference on Gravitational Waves, astro- 
ph/9911532| (1999). 



C. L. Fryer, W. Benz, S. A. Colgate and M. Herant, As- 
trophys. J., 516, 892 (1999). 

L. Rezzolla, F. K. Lamb, D. Markovic and S. L. Shapiro, 
Phys. Rev. D, 64, 104013 (2001). 

T. W. Baumgarte, S. L. Shapiro and M. Shibata, Astro- 
phys. J. Lett., 528, 29 (2000). 
S. L. Shapiro, Astrophys. J., 544, 397 (2000). 



17 



